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Abstract. We study multiple period states of a two-component unpolarized 
superfluid Fermi gas in an optical lattice along the Bardeen-Cooper-Schrieffer (BCS) 
to Bose-Einstein condensate (BEC) crossover. The existence of states whose period is 
a multiple of the lattice spacing is a direct consequence of the non-linear behavior 
of the gas, which is due to the presence of the order parameter associated with 
superfluidity. By solving Bogoliubov-de Gennes equations for a superfluid flow with 
finite quasimomentum, we find that, in the BCS side of the crossover, the multiple 
period states can be energetically favorable compared to the normal Bloch states and 
their survival time against dynamical instability drastically increases, suggesting that 
these states can be accessible in current experiments, in sharp contrast to the situation 
in BECs. 


PACS numbers: 03.75.Ss, 67.85.De, 03.75.Lm, 67.85.Hj 
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1. Introduction 

Density structures and patterns caused by the interplay of competing effects are 
ubiquitous in nature. Examples are the competition between dispersion and non¬ 
linearity, which yields solitons jT], the competition between crystalline order and Peierls 
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instability of conduction electrons, which results in charge and spin density waves [2], 
or the Fulde-Ferrcll-Larkin-Ovchinnikov (FFLO) state [3], with a spatially-dependent 
pairing held originating from the competition between the mismatch of the Fermi 
surfaces in the imbalanced systems and the energy gain by the condensation; but similar 
conditions also occur in the “pasta” phases in neutron stars [3], in nuclear halo [5], and 
in superfluid 4 He [6]. In the case of superhuids in a periodic lattice, non-linearity due 
to the presence of the order parameter, which favors a quadratic energy dispersion, can 
lead to the persistence of the quadratic-like dispersion beyond the Brillouin zone edge 
and give rise to non-trivial loop structure called “swallowtail” in the energy band [THUj. 
Due to their high controllability jT5lfl6] . ultracold gases offer an excellent test bed for 
exploring these intriguing phenomena. 

For atomic Bose-Einstein condensates (BECs) flowing in periodic potentials with 
finite quasimomentum, it was found that non-linearity of the interaction can cause the 
appearance of stationary states whose period is not equal to the lattice constant as in 
the usual Bloch states, but is a multiple of it [T6!fT8] : such states are called multiple (or 
n-tuple) period states. In BECs without long-range interaction, however, these states 
are energetically unfavorable compared to the normal Bloch states and unstable against 
small perturbations na. Here we investigate multiple period states in atomic Fermi 
superfluids, which are particularly interesting for their analogs in condensed matter 
physics and nuclear physics, such as superconducting electrons in solids and superfluid 
neutrons in neutron stars jl9[|20]. Furthermore, by using Feshbach resonances one can 
continuously go from the Bardeen-Cooper-Schrieffer (BCS) to the BEC regimes |21.22!. 
thus allowing one to understand Bose and Fermi superfluids from a unified perspective. 
Unlike the case of Bose gases, little has been studied about multiple period states in 
Fermi gases and their existence itself along the BCS-BEC crossover is an open question. 
In this work, we show that they indeed exist and they can be energetically favorable 
compared to the normal Bloch states in the BCS regime. Furthermore, we find that, 
despite being dynamically unstable, their lifetime becomes drastically long by going 
toward the deep BCS limit, possibly allowing for their experimental observation. 

This paper is organized as follows. After we explain the basic formalism employed 
in the present work in section [2l we show that multiple period states appear in Fermi 
superfluids along the BCS-BEC crossover and discuss their stationary properties in 
section [3] and [4j We then discuss their dynamical stability in section [5j Finally, this 
paper is concluded in section [6] 

2. Setup and basic formalism 

We consider an equally populated (unpolarized) two-component Fermi gas in the 
superfluid phase at zero temperature, moving in a one-dimensional (ID) optical lattice, 

I4ct(r) = V ext (z) = V 0 sin 2 q B z = sE R sin 2 q B z , (1) 

where s is the dimensionless parameter of the lattice height, E R = h 2 q 2 B /2m is the recoil 
energy, m is the atom mass, q B = Tr/d is the Bragg wave vector, and d is the lattice 
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constant. Note that qs differs from the fundamental vector of a ID reciprocal lattice, 
2i\ jd, by a factor of 2. The gas is uniform in the transverse directions and we look for 
stationary states of the system in the BCS-BEC crossover by numerically solving the 
Bogoliubov-de Gennes (BdG) equations [221123] : 

H{ r) A(r) \ / ufl r) 

A*(r) -H( r) J vflr) 



where H( r) = —h 2 V 2 /2m + 14 xt (r) —/i, ufl r) and 7y(r) are quasiparticlc amplitudes, and 
e % is the corresponding quasiparticle energy. The chemical potential p is determined from 
the constraint on the average density no = N/V = V^ 1 f n(r) dr = 2V~ 1 Yfl f |vj(r)| 2 dr 
with N being the number of particles and V being the volume, and the pairing field 
A(r) should satisfy a self-consistency condition A(r) = — g JA uflr)v*(r), where g is the 
coupling constant for the s-wave contact interaction which needs to be renormalized 
j24|f27] . The total energy E is given by 


E 



+ Vext(r)n(r) + —| A(r) | 2 
9 


In this formalism, a stationary motion of the superfluid in the ^-direction, relative 
to the infinite periodic potential at rest, is described by solutions of equation (J2]) 
with quasimomentum P per atom (not per pair), or the corresponding wave vector 
Q = P/h , such that the quasiparticle amplitudes can be written in the Bloch form 
as ufl r) = uflz)e l ® z e l]i ' r and vflr) = v t (z)e~ l ® z e tk ' r leading to the pairing field as 
A(r) = e t2 ® z A(z). Here A(z), uflz), and vflz) are complex functions with period u 
times d, with v G {1,2,3, • • •}, and the wave vector k z lies in the first Brillouin zone 
for a superccll (a cell containing several primitive cells) with period v, i.e., \k z \ < qslv. 
This Bloch decomposition transforms (EJ into the following BdG equations for uflz) and 
vflz): 

H q (z) A (z) \ / uflz) 

A *{z) -H-q(z) ) y vflz) 

where 

Hq(z) = \k'_]_ + (—id z + Q + k z ) 2 ] + Kxt (z) — p . 

Here, k\ = k 2 -\-k 2 and the label i represents the wave vector k as well as the band index. 
We solve this BdG equations (J3]) for a supercell with period v (—vd/2 < z < ud/2) to 
obtain the period-z/ states. As in P2112D1EB1E9] (see also appendix iv), the detailed 
procedure follows these steps: Starting from an initial guess of A (z) and p (final results 
are robust to the choice of the initial guess), we diagonalize the matrix of the left-hand 
side of the BdG equations fl3]) and obtain e il hi, and ig. Based on the obtained and 
we calculate the average number density n 0 and A(z). If the resulting n 0 does not 
agree with a given target value n{ argct , we update p according to the difference between 
these values. Specifically, in our calculations, we set the updated value p' using the 
following formula: p! — p ( n ) ) arget / no ) ?) with rj < 1 such as 77 = 2/3, 1/3, 1/5, etc. 
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Figure 1 . Profiles of (a) the magnitude of the pairing field |A(z)| and (b) the 
density n(z) of the lowest period-doubled states with period 2d (period-2 states) at 
their Brillouin zone edge P = P e d ge /2 = hqs/4: for three different values of l/kpa s : 
1/kpCLs = —1 (red solid line), —0.5 (green dashed line), and 0 (blue dotted line). They 
are obtained for s = 1 and Ef/Er = 0.25. At P = P e d ge /2, A(z) of the period-2 
states has a node. 


Until fi converges within 10 ~ 7 E Rl we iterate the above procedure using the obtained 
A(z) and updated /i. We check that other key quantities also converge with sufficient 
accuracy. 

In the following, we mainly present the results for s = 1 and 2 with Ep/Ep = 0.25 
as examples, where Ep = K 2 k 2 F /{2rn) and kp = are the Fermi energy and 

wavenumber of a uniform free Fermi gas of density n 0 . These values fall in the range 
of parameters of feasible experiments We have performed systematic calculations 
for different values of Ep/Ep (0.2, 0.25, and 0.5) as well and checked that the main 
results of stationary properties remain qualitatively the same. We denote by P e d ge the 
quasimomentum P per atom at the edge of the Brillouin zone for the normal Bloch 
states with period 1. For superfluids of fermionic atoms, P e d g e = hq B / 2. Note that this 
value differs by a factor of 2 from that of superfluids of bosonic atoms, hq B , because the 
elementary constituents of Fermi superfluids are pairs of fermionic atoms. 

3. Stationary solutions 

In figures [0(a) and ED)b) , we show the profiles of the pairing field |A(z)| and the number 
density n(z) of the lowest period-doubled states (period-2 states), respectively. Here we 
set P = P e dge/2 = hq B / 4 at the Brillouin zone edge of the period-2 states, where the 
feature of the period-2 states appears most prominently 0. The feature of the period 
doubling and the difference between the regions of— 1 < z/d < 0 and 0 < z/d < 1 can 

f For the BEC case it has been shown that, with increasing non-linearity (i.e., the interaction strength 
gk) from the linear limit (gi, = 0), the period-doubled states start to appear at the Brillouin zone edge 
of the period-2 system and their band extends in the Brillouin zone El. In this sense, the Brillouin 
zone edge for the period-2 system is a representative point for period-doubled states. 
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Figure 2. Profiles of (a) the pairing field | A(z)| and (b) the density n(z) of the normal 
Bloch states (magenta dashed line for l/k F a s = —1 and cyan dashed-dotted line for 
l/k F a s = 0) and period-doubled states (red solid line for 1 /kpa s = —1 and blue dotted 
line for 1 /kpa s = 0) at P = P e< j ge /2. Here we set s = 1 and E F /E R = 0.25 as in 
figure Q] 


be clearly seen in |A(z)| at any value of l/k F a s . At P = P e d ge /2, A(z) of the period- 
2 states has a node [see z/d = 0.5 in figure [0(a)] and consequently the supercurrent 
is zero ( dpE = 0; see figure [3]). On the other hand, especially in the deeper BCS 
side (1 /k F a s = —1), the difference in n(z) between the regions of —1 < z/d < 0 and 
0 < z/d < 1 is small [see the red line in figure H^b)]. Around z/d = 0.5, where |A(z)| 
vanishes, the density remains large, suggesting the existence of Andreev-like localized 
states. The density difference between the two regions becomes larger with increasing 
1 /k F a s toward the BEC regime. 

In figure El we compare the profiles of |A(^)| and n(z) between the normal Bloch 
states and the period-doubled states for 1 /kpa s = —1 and 0. Since we set P at the 
Brillouin zone edge of the period-doubled states (P = P e d ge /2), where their supercurrent 
is zero, A(z)’s of the period-doubled states have a node while those of the normal Bloch 
states do not. Note that, in the BCS regime of 1 jkpa s = —1, n(z)’s of the normal Bloch 
state and the period-doubled state are almost the same, but |A(z)|’s are significantly 
different. 

It is instructive to consider the deep BCS limit of l/k F a s —* — oo. There, A(z), 
which is the origin of the non-linearity, vanishes and n(z) of the neighboring sites 
becomes identical so that the nature of the period doubling disappears in this limit. 
We observe that by going to the deep BCS regime, where A(z) and the supercurrent 
are infinitesimally small, the energy difference E(P) —E( 0) for period-2 states decreases 
[i.e., E(P) becomes more flat] and period-2 states at P = P e d ge /2 approach the 
normal Bloch state at P = 0. These observations are consistent with the fact that, 
if A(z) = 0, our non-linear BdG equations reduce to the linear Schrodinger equation, 
whose solutions have the periodicity of the lattice due to the Bloch theorem. Multiple 
period states are hence possible only in the superfluid phase. It is worth mentioning 
here that these multiple-period states are essentially different from the FFLO [3] or 
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Figure 3. Energy E per particle in units of Er as a function of the quasimomentum 
P. Here we set s = 1, Er/Er = 0.25, and 1 /kpa s = —1. The normal Bloch states 
with period d are shown by the blue dotted line with • symbols, and the multiple 
period states are shown by the red solid line with + (period 2), the green dashed line 
with O (period 3), the purple dashed-dotted line with O (period 4), and the magenta 
solid line with □ (period 5). The dotted lines at P/P e d ge = 1/3, 0.25, and 0.2 show 
the first Brillouin zone edge for the period-3, 4, and 5 states, respectively. The inset 
shows that the lowest band of the period-5 states continues beyond the first Brillouin 
zone edge and the swallowtail appears (shown by the red arrow). 


soliton lattice [31] states in the imbalanced (polarized) systems. In onr case, the non¬ 
trivial spatial dependence of the pairing field is a purely non-linear phenomenon caused 
by the presence of the superfluid order parameter, while in the other cases it is due 
to the non-zero center-of-mass momentum of the pair, which requires the mismatch of 
the Fermi surfaces between two components. The multiple period states studied in the 
present work is also different from the charge density wave due to the nesting of the 
Fermi surface (see appendix i for details). 

As a final comment on the spatial structure of stationary solutions, it is worth 
noting that the presence of a node in A is a sufficient condition for a zero supercurrent, 
but it is not a necessary condition. For example, at the Brillouin zone center, the 
supercurrent is of course zero because the phase of A is constant (P = 0), but A does 
not have a node. On the other hand, for nonzero P, the phase of A depends on the 
position. Therefore, when the supercurrent is zero at the Brillouin zone boundary, A 
must have a node. States with more nodes in A have higher energy in general. As to 
the lowest periodic states which we discuss in the present work, the number of nodes 
is thus one per superccll at the Brillouin zone boundary. To minimize the energy, the 
node is located at the potential maximum. 
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4. Energetics 

In figure [3l we plot the energy per particle of the lowest band as a function of the 
quasimomentum P for the normal Bloch states (blue dotted line) and the multiple 
period states. We show the results at 1 /kpa s = —1 in the BCS side. In the region of 
small P, all lines of the multiple period states collapse onto the line of the normal Bloch 
states, as they are all equivalent in this region, the states with period 1 being just a 
subset of all the multiple period states [§[ 

Conversely, the multiple period states for small v (z/ < 4 in the case of figure [3]) 
become energetically more stable than the normal Bloch states near the first Brillouin 
zone edge of each multiple period state, i.e., P < P e d ge /z/ for period-z/ states. In 
particular, the period-doubled states are the lowest in energy in a wide range of P 
[In the case of figure [3l period-2 states are always energetically lower than period-3 
states, which holds in the region of 0.245 < Ep/Ep < 0.4 for the same values of s = 1 
and 1 /kpa s = —1 (see appendix ii for details).]. This is in striking contrast to the 
situation in BECs and in the BEC regime of the BCS-BEC crossover (see later), where 
the lowest band of normal Bloch states is always lower in energy than the multiple period 
states [Eli the latter appear as an upper branch of the swallowtail band structure (with 
d 2 P E > 0) around the Brillouin zone edge of the respective multiple period states m. 
Figure [3] shows instead that, in the BCS regime, the lowest band of the multiple period 
states for small v has dpE < 0 near the Brillouin zone edge P < P e dge/^- 

At first sight, this seems to imply a pathological situation in which //-period states 
with large v would be lower in energy than the normal Bloch states even in the limit of 
P —)■ 0. However, this pathological situation is saved by the emergence of the swallowtail: 
The multiple period states with large v continue being almost identical to the normal 
Bloch states, and keep their nearly quadratic dispersion around P ~ 0 even beyond 
their first Brillouin zone edge at P e dge/zfi which results in the swallowtail band structure 
for the period-z/ system. In the case of figure [31 the swallowtail starts to appear at v = 5 
(see the inset of figure [3]) . 

In figure [4J, we show the total energy difference A E = P 2 — Pi between the normal 
Bloch states (Pi) and the period-2 states (P 2 ) at P = P e dge/2 = hq B / 4 along the BCS- 
BEC crossover. As we have seen in figure El the period-doubled states are energetically 
more stable (i.e., A E < 0) in the BCS regime. With l/kpa s increasing from the deep 
BCS limit, AP increases from a negative value and finally period-doubled states become 
higher in energy than normal Bloch states (i.e., AP > 0) in the BEC side. Here we 
point out that, in the region of AP < o, the period-doubled states form a band which is 
convex upward and smoothly connects to that of the normal Bloch states (see figure [3]) , 
and the swallowtail does not exist. On the other hand, in the region of AP > 0, they 

§ The critical quasimomentum P c for the pair-breaking instability of the normal Bloch state is 
0.147P e dge in this case. Note that P c seems to coincide with the value at which the multiple period 
states start to separate from the normal Bloch state. This suggests that emergence of the multiple 
period states lead to the Landau instability of the normal Bloch state. 
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1 / k F a s 

Figure 4. Difference A E = E 2 — E\ of the total energy per particle in units of Er 
between the normal Bloch states {Ef) and period-doubled states (E 2 ) at P = P e d ge /2. 
The parameters we have used are s = 1, 2 and Er/Er = 0.25. The red solid line with 
+ (s = 1) and blue solid line with x (s = 2) show the results obtained by solving the 
BdG equations and the green dashed line shows the results by the Gross-Pitaevskii 
equations for parameters corresponding to s = 1 and Ef/Er = 0.25. 


form a concave upper edge of the swallowtail, which is located above the crossing point 
(“x’’-like shape) of the swallowtail. Therefore, hysteresis caused by the swallowtail, 
which could be observed in the latter region (A E > 0) of the BCS-BEC crossover and 
in BECs [ll|, would disappear in the former region (A E < 0). 

We also show the results in the BEC side obtained by solving the Gross-Pitaevskii 
(GP) equation for corresponding parameters (the green dashed line). Namely, GP 
equation for bosons of mass mb = 2m interacting with scattering length = aa s = 2 a s 
for the mean-held theory in an optical lattice 2V ex t(z). We can relate l/kpci s and 
gbnb/E Rb , where gb = ArafFablmb, nb = nj2 , and E Rb = h 2 (n/d) 2 /2mb are the 
interaction parameter, the average density, and the recoil energy of bosons, as l/k R a s = 
(4a/37r) (gbUb/ -Er^) -1 (Er/Er). We note that A E of the GP results approaches zero 
from positive values with increasing 1 /k F a s . This suggests that, before the BdG results 
(the red and blue solid lines) converge to the GP ones in the deep BEC regime, A E 
takes a maximum value. 

For different strength s of the lattice, we see that the curve of A E is somewhat 
shifted towards the BCS side with increasing s, so the period-doubled states become 
less stable (see the red solid line with + for s = 1 and the blue solid line with x for 
s = 2 in figure 0]) . This might be due to the formation of bosonic molecules of fermionic 
atoms induced by the external lattice potential [28, 32, 33]. 

The energetic stability of multiple period states in the BCS regime can be physically 
understood as follows. Let us consider the different behavior of A(z) and n(z) for a 
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period-2 state and a normal Bloch state at P = P e d ge /2. In the case of a normal 
Bloch state, since |A(z)| is exponentially small in the BCS regime, we can distort 
the order parameter A( 2 ) to produce a node, like the one in the period-2 state, with 
a small energy cost (per particle) up to the condensation energy \E cond \/N <C E F , 
where E cond = g~ l f d 3 r |A(r)| 2 . However, making a node in A(^) kills the supercurrent 
j = V~ 1 dpE, which yields a large gain of kinetic energy (per particle) of the superfluicl 
flow of order rs_/ P edge/ m E r . Even if A( 2 ) is distorted substantially to have a node, 
the original density distribution of the normal Bloch state is almost intact so that the 
increase of the kinetic energy and the potential energy clue to the density variation is 
small. Therefore, the period-2 state is energetically more stable than the normal Bloch 
state in the BCS regime. In the above discussion, the key point is that A(z) and n(z) can 
behave in a different way in the BCS regime. On the other hand, in the BEC limit, the 
density is directly connected to the order parameter as n(z) oc |A(z)| 2 , and distorting 
the order parameter accompanies an increase of the kinetic and potential energies clue 
to a large density variation. 

More generally, for period-// states in comparison with a normal Bloch state at 
P = -P e dge/O the energy cost to distort A(z) to have a node is up to ~ \E cond \/N <C Ep, 
but the energy gain is of order rs./ P edge / m " 2 E R /u 2 , which is reduced by a factor of 
v 2 . We thus see that period-// states with sufficiently large v cannot be energetically 
more stable than the normal Bloch state as has been observed before. 

5. Dynamical stability and survival time 

So far, we have seen that multiple-period states exist as energetically stable stationary 
solutions of the BdG equations. The next important issue is their dynamical stability, 
that is, whether and how long they can survive under small perturbations, which are 
unavoidable in experiments. We face this problem by performing numerical simulations 
based on the time-dependent BdG (TD-BclG) equations. 

A crucial difference between the stationary calculations (time-independent BdG) 
in the previous sections and the dynamical (time-dependent BdG) calculations in this 
section is the following. The ideal configuration to study the stationary solutions with a 
given periodicity v is a supercell with v sites under the Bloch-wave boundary conditions, 
as we have done in the previous sections. Conversely, dynamical calculations has to 
account for excited states with any wavelength, possibly including long wavelength 
perturbations which may trigger a dynamical instability, so that the Bloch-wave 
boundary conditions cannot be used. We instead solve the TD-BdG equations in a large 
computational box of length L z in the ^-direction, including a sufficiently large number 
of supercells of v sites, with the periodic boundary conditions, in order to mimic an 
infinite system with good enough accuracy. We use L z = 32 d to 64 d] L z is chosen to be a 
multiple of 8 d for convenience, so that the allowed values of the wave vector k z discretized 
as A k z = 27t/L 2 are commensurate with the value of the quasimomentum at the first 
Brillouin zone edge for the period-2 states, P edge /{2n) = q R /4. Finally, stationary BdG 
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Figure 5. Time evolution of the magnitude of the pairing field |A(z, i)|/| A(z = 0, t = 
0)| of the period-doubled state at P = P e d ge /2 for 1/fc^a s = —1 (left panel), 0 (middle 
panel), and 0.5 (right panel). Here, s = 1 and Ef/Er = 0.25. Actual calculation has 
been done for L z = 32 d in the cases of 1 /fc_pa s = 0.5 and 0 and for L z = 48d in the 
case of 1/kpcis = —1; a part of the system is shown in the figure. 


calculations for such large boxes with the periodic boundary conditions are not feasible 
with our current computational resources, since they require long iterative procedures 
for many values of P; a direct comparison between stationary and time-dependent BdG 
results would be possible only for smaller values of L z , corresponding to less than about 
ten lattice sites, for which the extrapolation to an infinite system would be unreliable. 

As initial configuration of the TD-BdG simulations, we use a configuration based on 
the stationary solution of the BdG equations, which is constructed as follows. Among the 
quasiparticle amplitudes rq and ig obtained by solving equation (J3]), we select those with 
(quasi)wave vectors k z equal to multiples of A k z = 2i t/L z . In this way, we construct the 
approximate stationary solution of the BdG equations (J2J) with the periodic boundary 
condition. Then we integrate the TD-BdG equations using a 4-tli order predictor- 
corrector method. The basic structure of the code is the same as the one in [34j. 

The right and middle panels of figure [5] show the time evolution of |A(z)| of the 
period-doubled state at P = P e d ge /2 in the BEG side (1 /kpa s = 0.5) and at unitarity, 
respectively. We see that |A(z)| [and n(z)\ does not keep its initial profile and large- 
amplitude oscillations triggered by the dynamical instability set on at t ~ 55 H/Er in 
the right panel and t — 130 h/E r in the middle panel. We also notice that the TD-BdG 
simulations allow us to identify the spontaneously growing excitations which trigger 
the instability. The wavelength of the growing mode is 4 d, 4 d, and 2d in the case of 
1 /kFa s = —1 (left panel of figure EJ) , 0 (middle panel), and 0.5 (right panel), respectively. 

It is remarkable that the survival time r surv of the period-doubled states until they 
are destroyed by the large-amplitude oscillations drastically increases as going toward 
the BCS side. In the left panel of figure [5l we show the time evolution of |A(z)| of 
the period-doubled state at P = P e d ge /2 for 1 /kpa s = —1. In this realization, the 
period-doubled state almost keeps its initial profile of |A(z)| until t — 900 h/E Rl and 
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Figure 6. Growth rate 7 of the fastest growing mode (black solid line) and survival 
time T surv of the period-doubled state at P = P e d ge /2 (s = 1 and E F /E R = 0.25). 
Blue dashed-dotted, green dotted, magenta dashed double-dotted, and red dashed lines 
show r surv for relative amplitude fj( 0 ) of the initial perturbation of 10 %, 1 %, 0 . 1 %, and 
0.01%, respectively. (The regions of l/k F a s < —2 and > 0 are extrapolated.) 


even longer for n(z) because only a small fraction of particles participate in the pairing 
in the BCS regime. 

To further analyze the time scale of the deviation \A(z, f)| — |Ao(z)| from the true 
stationary state A 0 (z), we take its spatial Fourier transform and look for modes with 
exponentially growing amplitudes |? 7 (i)| = |r/(0)| e 7 *. From a fit we extract the growth 
rate 7 of the fastest growing mode. The growth rate 7 corresponds to the imaginary part 
of the complex eigenvalue for the fastest growing mode obtained by the linear stability 
analysis [[I61I35J. This is intrinsic property of the initial stationary state independent of 
the magnitude of the perturbation. 

The resulting 7 is shown by the black solid line in figure [ 6 ] which clearly shows 
the suppression of 7 with decreasing l/k F a s . In practice, the survival time r surv of the 
period-doubled states depends on the accuracy of their initial preparation. We estimate 
r surv with ?)( 0 )e 7t rs j 1 , where 7 ( 0 ) is the relative amplitude of the initial perturbation 
with respect to |A 0 | for the fastest growing mode. In figure El we show r surv for four 
values of 77 (C)). This result suggests that if the initial stationary state is prepared within 
an accuracy of 10 % or smaller, this state safely sustains for time scales of the order of 
100 Ti/Er or more in the BCS side, corresponding to r surv of more than the order of a few 
milliseconds for typical experimental parameters [3D]: For E fi h = 2n x 7.3kHz x h used in 
the experiment of [3D], 1 h/E R = 0.0109 msec. In the deep BCS regime (1 /k F a s -C — 1), 
T aurv increases further and may become larger than the time scale of the experiments, so 
that the period-doubled states can be regarded as long-lived states and, in addition, since 
they have lower energy than the usual Bloch states in a finite range of quasimomenta, 
they could be realized by, e.g., quasi-adiabatically increasing P from the ground state 
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at P = 0, which is the normal Bloch state. 

Finally, it is worth noting that the BCS transition temperature T c is roughly 
estimated as T c T F e n ^ 2kFas with T f = Ep/ks and k F is the Boltzmann constant. 
For the above value of E n used in the experiment of [30], T c ~ 200nK at 1 /k F a s = — 1, 
50nK at 1 /k F a s = —2, and lOnK at l/k F a s = —3. Therefore, superfluidity can be 
realized in the whole region shown in figure [6] in the current experiments. 

6. Conclusion 

We have studied multiple period states, especially period-doubled states, of superfluid 
Fermi gases in an optical lattice and have found that they can be energetically more 
stable than the normal Bloch states and their survival time can be drastically enhanced 
in the BCS side. The multiple period nature distinctly appears in the pairing field, 
which could be observed by the fast magnetic field sweep technique [361137] , It is also 
interesting to point out that the emergence of the period-doubled states in the BCS side 
is closely connected to the disappearance of the swallowtails, which exist in the BEC 
side (section |4j) . As a consequence, hysteresis of the superfluid circuits [ITj , which could 
be observed experimentally in the BEC side, would disappear by sweeping to the BCS 
regime. We hope our work will stimulate future experimental studies. 
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Appendix 

In this appendix we provide additional information about the momentum distribution, 
comparison of the energy between period-2 and -3 states, comments on rational and 
irrational number periods, and detailed information of the numerical calculations. 

i) Momentum distribution 

The period n-tupling studied in the present work is different from the charge density 
wave due to the nesting of the Fermi surface. In figure [7] we show the momentum 
distribution n(k) of the normal Bloch states at P = 0. At l/k F a s = —1 [panel (a)], 
there is a plateau around k± = k z = 0 with a smeared Fermi surface whose width 
is characterized by ~ |A| 1//2 , while, at l/k F a s = 0 [panel (b)], n(k) shows a peak at 
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(b) 



kj_/k F 


Figure 7. Momentum distribution n(k)kp at P = 0 in the BCS side of 1 /k F a s = — 1 
[panel (a)] and at unitarity 1 /k-pa s = 0 [panel (b)] for s = 1 and Ep/Ep = 0.25. The 
vertical axis shows the wave vector k z in z-direction and the horizontal axis shows 
the wave vector kj_ in the transverse directions. Contours at n(k)kp = 0.2 (green 
dashed-dotted), 0.4 (blue dashed), 0.6 (red solid), and 0.8 (cyan dashed) are shown. 



Figure 8. Same as figure [7] for a stronger periodic potential with s = 2. As in 
figure[7|a), we set P = 0, \/k-pa s = —1, and Ep/Ep = 0.25. 


k± = k z = 0 rather than a plateau. Note that even though the system is in a periodic 
potential with nonzero s, the Fermi surface is almost spherical. 
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Figured] is the same as figure [3(a), but for a stronger periodic potential with s = 2. 
The plateau region of the momentum distribution n(k) is significantly smaller compared 
to that for s = 1 shown in figure [7j(a) . This is due to the formation of bound bosonic 
dimers induced by the stronger periodic potential [281 1321133] . However, also in this case, 
n(k) is almost isotropic although it is more compressed in the fc_i_-directions compared 
to the case of s = 1 [figure [3(a)] . 

ii) Comparison of the energy between period-2 and -3 states 



Figure 9. Energy difference E 2 — E 3 between the period-2 and -3 states at 
P = P e dge/3, where E v (u = 2,3) represents the energy of the period-^ state. Here we 
set s = 1 and 1 /fc_pa s = — 1, which are the same as in figure [3] 

In figure [9j we show the energy difference E 2 — E 3 between the period-2 and -3 
states at the Brillouin zone edge of the period-3 states, P = P e d g e/3. For the parameter 
of figure [3] Ef/Er = 0.25, we see that E 2 — E 3 < 0 and thus the period-2 states are 
energetically lower than the period-3 states in the whole Brillouin zone of the latter (see 
figure [3]) . This figure shows that the period-2 states are always energetically lower than 
the period-3 states in the region of 0.245 < Er/Er < 0.4 for s = 1 and 1 /kpa s = —1. 

iii) Comments about rational and irrational number periods 

In principle, multiple period states with rational number periods are covered by 
our calculations with a supercell. Specifically, using a supercell with period is, we can 
describe period-(is/rj) states, where v and g are natural numbers. On the other hand, 
states with irrational number periods are excluded, which are beyond the scope of 
the present study. In our numerical calculations, which cover multiple period states 
with rational number periods, neither states with is < r] nor states whose period is 
incommensurate with the lattice constant appear as the lowest energy state. Therefore, 
it is probable that the multiple period states with irrational number periods could not 
be the energetically minimum states. 
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iv) Detailed information of the numerical calculations 

We set the parameters for the numerical calculations such as the number of grid 
points depending on the system parameter values to ensure the convergence. Here 
we provide detailed parameter values for the numerical calculations for s = 1 and 
Ef/Er = 0.25 as a typical example. 

In the transverse directions, we impose periodic boundary conditions with a large 
box size with L F /d = 24. Regarding the calculations of the stationary BdG equations, 
the cutoff energy E c we use is, for example, E c /Ep = 40 for 1 /k F a s = —1, 40-200 
(mainly 100) for 1 /kpCh s = 0, and 100-200 for l/k F a s = 0.5. The number of the 
grid points for k z within the first Brillouin zone — qs/v < k z < qs/v for period v 
is 100, 100, 75, 50, and 40 for v — 1, 2, 3, 4, and 5, respectively. The number 
of the grid points in a supercell in the 2 direction, —v/2 < z/d < z^/2, is 200 
for v = 1 and lOOzv for the other values of v. Regarding the TDBdG simulations 
shown in figure [5j the time discretization At is 0.00008/1/77}? = 0.00032 h/E R for 
l/k F a s = — 1 (left panel), 0.00005/1/77}? = 0.0002 H/Er for 1 /k F a s = 0 (middle panel), 
and 0.00004 h/E F = 0.00016/1/77/} for l/k F a s = 0.5 (right panel). Discretization of 
k z and z are A k z /qB = 0.025 and Az/d = 0.025, respectively. Throughout the time 
evolution of figure [5J the total number of particles is conserved perfectly within the 
significant digits and the energy is conserved within 0.045% for l/k F a s = —1 (left 
panel), 0.74% for 1 /k F a s = 0 (middle panel), and 2.2% for 1 /k F a s = 0.5 (right panel). 
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